Everything is related to everything else, but near things are more related than distant things.
- Waldo Tobler
First Law of Geography
This to be wary of
- Does this actually need to be on a map?
This to be wary of
- Does this actually need to be on a map?
- Is the structure in my map just another variable lurking in the background?
Recreating Snow's Map
Step 1
Pull in cholera deaths and pump location data.
Notes:
- Will involve shapefiles, which can be points, paths, lines (vector), or polygons.
- Based on ESRI's formatting (they develop ArcView GIS)
- Uses
spandrgdalpackages.
DEMO
## OGR data source with driver: ESRI Shapefile ## Source: "/Users/andrewbray/Dropbox/Teaching/2018-spring/math-241/course-materials/data/snow/", layer: "Cholera_Deaths" ## with 250 features ## It has 2 fields
## OGR data source with driver: ESRI Shapefile ## Source: "/Users/andrewbray/Dropbox/Teaching/2018-spring/math-241/course-materials/data/snow/", layer: "Pumps" ## with 8 features ## It has 1 fields
Cholera data
cholera_locs <- coordinates(cholera) head(cholera_locs)
## coords.x1 coords.x2 ## [1,] 529308.7 181031.4 ## [2,] 529312.2 181025.2 ## [3,] 529314.4 181020.3 ## [4,] 529317.4 181014.3 ## [5,] 529320.7 181007.9 ## [6,] 529336.7 181006.0
class(cholera_locs)
## [1] "matrix"
cholera_locs <- as.data.frame(cholera_locs)
library(tidyverse) ggplot(cholera_locs, aes(x = coords.x1, y = coords.x2)) + geom_point()
Recreating Snow's Map
Step 2
Pull in map of London. Superimpose cholera and pump data.
Notes:
- A map involves data originally collected on the sphere.
- We need both the basemap as well as the projection that will determine how it is displayed on our 2D screen.
- We'll use
ggmapbut there are several other packages to do this step.
DEMO
London basemap
library(ggmap)
m <- get_map("John Snow, London, England", zoom = 17, maptype = "roadmap")
ggmap(m)
ggmap(m) + geom_point(cholera_locs, aes(x = coords.x1, coords.x2))
This won't work: units down align.
head(cholera_locs, n = 2)
## coords.x1 coords.x2 ## 1 529308.7 181031.4 ## 2 529312.2 181025.2
str(m)
## 'ggmap' chr [1:1280, 1:1280] "#FAF5EA" "#FDF9F2" "#FDF9F2" "#FDF9F2" ... ## - attr(*, "bb")='data.frame': 1 obs. of 4 variables: ## ..$ ll.lat: num 51.5 ## ..$ ll.lon: num -0.14 ## ..$ ur.lat: num 51.5 ## ..$ ur.lon: num -0.133 ## - attr(*, "source")= chr "google" ## - attr(*, "maptype")= chr "roadmap" ## - attr(*, "zoom")= num 17
Projections
A projection transforms coordinates in 3D space to coordinates in 2D space.
library(maps)
map("world", project = "mercator", wrap = TRUE)
## Warning in map("world", project = "mercator", wrap = TRUE): projection
## failed for some data
map("world", projection = "cylequalarea", param = 45, wrap = TRUE)
Projection types
In general, projections can focus on preserving one of the following (at the expense of the other).
- Shape/Angle: Mercator, Lambert
- Area: Gall-Peters, Albers
Projection information and other info about the geographic data is monitored via a standardized coordinate reference system (CRS) using a library called PROJ.4.
CRS
library(sp) proj4string(cholera)
## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
Standard code formats:
- EPSG:4326 (WGS84) used by GPS and Google Earth
- ESPG:3857 Used for map tiles in Googe Maps, Open Street Maps, etc.
- ESPG:27700 (OSGB 1936) commonly used in Britain
CRS
The basemap is actually in EPSG:4326, so let's project our cholera data into those coordinates.
proj4string(cholera)
## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
cholera_proj <- cholera %>%
spTransform(CRS("+init=epsg:4326")) %>%
as.data.frame()
cholera_proj
## Id Count coords.x1 coords.x2 ## 1 0 3 -0.1363245 51.51291 ## 2 0 2 -0.1362775 51.51285 ## 3 0 1 -0.1362473 51.51281 ## 4 0 1 -0.1362063 51.51275 ## 5 0 4 -0.1361612 51.51269 ## 6 0 2 -0.1359313 51.51267 ## 7 0 2 -0.1365949 51.51285 ## 8 0 2 -0.1364394 51.51282 ## 9 0 3 -0.1366706 51.51281 ## 10 0 2 -0.1366179 51.51292 ## 11 0 2 -0.1367315 51.51287 ## 12 0 1 -0.1369572 51.51295 ## 13 0 3 -0.1368201 51.51271 ## 14 0 1 -0.1367729 51.51266 ## 15 0 4 -0.1367316 51.51261 ## 16 0 1 -0.1370399 51.51273 ## 17 0 1 -0.1370924 51.51265 ## 18 0 1 -0.1363187 51.51267 ## 19 0 4 -0.1362593 51.51260 ## 20 0 3 -0.1362059 51.51255 ## 21 0 2 -0.1371564 51.51293 ## 22 0 1 -0.1371933 51.51308 ## 23 0 2 -0.1374392 51.51289 ## 24 0 2 -0.1373643 51.51287 ## 25 0 2 -0.1372579 51.51290 ## 26 0 1 -0.1371470 51.51313 ## 27 0 1 -0.1372026 51.51318 ## 28 0 3 -0.1372505 51.51324 ## 29 0 1 -0.1372814 51.51317 ## 30 0 1 -0.1376339 51.51308 ## 31 0 1 -0.1377158 51.51315 ## 32 0 1 -0.1377106 51.51299 ## 33 0 2 -0.1380105 51.51307 ## 34 0 2 -0.1381139 51.51303 ## 35 0 1 -0.1384685 51.51279 ## 36 0 1 -0.1374886 51.51278 ## 37 0 1 -0.1380918 51.51250 ## 38 0 1 -0.1377217 51.51246 ## 39 0 2 -0.1377120 51.51238 ## 40 0 8 -0.1375814 51.51245 ## 41 0 2 -0.1374305 51.51252 ## 42 0 1 -0.1376037 51.51252 ## 43 0 1 -0.1368215 51.51232 ## 44 0 1 -0.1370186 51.51237 ## 45 0 1 -0.1364902 51.51202 ## 46 0 1 -0.1364290 51.51196 ## 47 0 4 -0.1363780 51.51192 ## 48 0 1 -0.1364593 51.51190 ## 49 0 1 -0.1365882 51.51200 ## 50 0 1 -0.1362121 51.51187 ## 51 0 1 -0.1360503 51.51194 ## 52 0 4 -0.1359787 51.51198 ## 53 0 1 -0.1360447 51.51186 ## 54 0 1 -0.1358441 51.51183 ## 55 0 1 -0.1357703 51.51185 ## 56 0 1 -0.1357216 51.51181 ## 57 0 1 -0.1353742 51.51203 ## 58 0 2 -0.1355748 51.51214 ## 59 0 1 -0.1354466 51.51218 ## 60 0 1 -0.1360897 51.51245 ## 61 0 1 -0.1359274 51.51226 ## 62 0 2 -0.1358134 51.51227 ## 63 0 1 -0.1357621 51.51222 ## 64 0 1 -0.1357193 51.51217 ## 65 0 2 -0.1359251 51.51240 ## 66 0 2 -0.1359561 51.51254 ## 67 0 1 -0.1358603 51.51256 ## 68 0 2 -0.1357804 51.51258 ## 69 0 3 -0.1357002 51.51261 ## 70 0 1 -0.1354832 51.51268 ## 71 0 4 -0.1353908 51.51270 ## 72 0 15 -0.1352533 51.51274 ## 73 0 3 -0.1351719 51.51276 ## 74 0 4 -0.1350995 51.51279 ## 75 0 5 -0.1348878 51.51265 ## 76 0 2 -0.1347238 51.51251 ## 77 0 1 -0.1348178 51.51241 ## 78 0 2 -0.1349176 51.51238 ## 79 0 1 -0.1349937 51.51235 ## 80 0 1 -0.1350932 51.51232 ## 81 0 1 -0.1352129 51.51227 ## 82 0 1 -0.1353670 51.51222 ## 83 0 1 -0.1347528 51.51236 ## 84 0 1 -0.1350242 51.51221 ## 85 0 1 -0.1349787 51.51214 ## 86 0 1 -0.1348168 51.51220 ## 87 0 1 -0.1347394 51.51211 ## 88 0 1 -0.1348317 51.51198 ## 89 0 1 -0.1347715 51.51194 ## 90 0 1 -0.1345911 51.51195 ## 91 0 1 -0.1345358 51.51190 ## 92 0 2 -0.1344959 51.51185 ## 93 0 1 -0.1344244 51.51176 ## 94 0 1 -0.1347040 51.51185 ## 95 0 4 -0.1343347 51.51148 ## 96 0 2 -0.1342521 51.51157 ## 97 0 1 -0.1341943 51.51152 ## 98 0 4 -0.1341110 51.51146 ## 99 0 4 -0.1335133 51.51137 ## 100 0 1 -0.1335386 51.51154 ## 101 0 4 -0.1337880 51.51174 ## 102 0 1 -0.1338034 51.51165 ## 103 0 1 -0.1338663 51.51170 ## 104 0 2 -0.1341597 51.51206 ## 105 0 1 -0.1342655 51.51207 ## 106 0 2 -0.1343699 51.51216 ## 107 0 3 -0.1344271 51.51222 ## 108 0 1 -0.1345094 51.51228 ## 109 0 4 -0.1345745 51.51234 ## 110 0 1 -0.1344772 51.51237 ## 111 0 1 -0.1345337 51.51243 ## 112 0 7 -0.1337228 51.51226 ## 113 0 3 -0.1335162 51.51233 ## 114 0 8 -0.1330388 51.51202 ## 115 0 1 -0.1329156 51.51169 ## 116 0 1 -0.1333610 51.51170 ## 117 0 5 -0.1334922 51.51264 ## 118 0 8 -0.1327876 51.51255 ## 119 0 2 -0.1328993 51.51266 ## 120 0 1 -0.1328308 51.51259 ## 121 0 1 -0.1329876 51.51273 ## 122 0 2 -0.1330345 51.51278 ## 123 0 1 -0.1331027 51.51287 ## 124 0 2 -0.1331502 51.51292 ## 125 0 2 -0.1336385 51.51297 ## 126 0 3 -0.1332914 51.51291 ## 127 0 1 -0.1335520 51.51302 ## 128 0 2 -0.1337383 51.51297 ## 129 0 2 -0.1334574 51.51308 ## 130 0 3 -0.1341951 51.51272 ## 131 0 1 -0.1341558 51.51267 ## 132 0 2 -0.1341344 51.51262 ## 133 0 1 -0.1340391 51.51254 ## 134 0 1 -0.1339959 51.51250 ## 135 0 1 -0.1338955 51.51237 ## 136 0 1 -0.1342263 51.51276 ## 137 0 3 -0.1344430 51.51295 ## 138 0 3 -0.1345338 51.51292 ## 139 0 3 -0.1346220 51.51289 ## 140 0 3 -0.1333934 51.51208 ## 141 0 1 -0.1331870 51.51207 ## 142 0 2 -0.1332902 51.51204 ## 143 0 1 -0.1333941 51.51201 ## 144 0 1 -0.1318768 51.51263 ## 145 0 1 -0.1316591 51.51272 ## 146 0 2 -0.1316900 51.51264 ## 147 0 1 -0.1313274 51.51275 ## 148 0 1 -0.1323920 51.51303 ## 149 0 1 -0.1324359 51.51312 ## 150 0 2 -0.1325501 51.51313 ## 151 0 2 -0.1324852 51.51301 ## 152 0 1 -0.1326663 51.51331 ## 153 0 1 -0.1326142 51.51321 ## 154 0 1 -0.1330981 51.51319 ## 155 0 1 -0.1331764 51.51332 ## 156 0 1 -0.1334045 51.51341 ## 157 0 1 -0.1333173 51.51309 ## 158 0 1 -0.1332787 51.51352 ## 159 0 5 -0.1326063 51.51338 ## 160 0 1 -0.1325292 51.51325 ## 161 0 1 -0.1327582 51.51356 ## 162 0 2 -0.1328407 51.51364 ## 163 0 2 -0.1328734 51.51369 ## 164 0 1 -0.1330518 51.51372 ## 165 0 1 -0.1327612 51.51381 ## 166 0 1 -0.1325726 51.51387 ## 167 0 1 -0.1325537 51.51385 ## 168 0 2 -0.1324627 51.51387 ## 169 0 2 -0.1324791 51.51389 ## 170 0 5 -0.1322153 51.51401 ## 171 0 1 -0.1323162 51.51399 ## 172 0 1 -0.1322437 51.51396 ## 173 0 1 -0.1321192 51.51399 ## 174 0 1 -0.1321392 51.51404 ## 175 0 4 -0.1320696 51.51405 ## 176 0 1 -0.1319574 51.51408 ## 177 0 2 -0.1318604 51.51407 ## 178 0 1 -0.1317867 51.51410 ## 179 0 1 -0.1328681 51.51533 ## 180 0 1 -0.1336530 51.51469 ## 181 0 1 -0.1337894 51.51464 ## 182 0 1 -0.1344164 51.51431 ## 183 0 3 -0.1351983 51.51433 ## 184 0 1 -0.1349775 51.51440 ## 185 0 1 -0.1340471 51.51399 ## 186 0 2 -0.1339720 51.51423 ## 187 0 1 -0.1332537 51.51396 ## 188 0 1 -0.1330836 51.51394 ## 189 0 1 -0.1332123 51.51434 ## 190 0 1 -0.1340977 51.51388 ## 191 0 1 -0.1339554 51.51389 ## 192 0 2 -0.1340428 51.51383 ## 193 0 1 -0.1338093 51.51371 ## 194 0 1 -0.1339705 51.51371 ## 195 0 2 -0.1337508 51.51364 ## 196 0 3 -0.1338695 51.51360 ## 197 0 2 -0.1346203 51.51385 ## 198 0 1 -0.1347223 51.51382 ## 199 0 1 -0.1346163 51.51403 ## 200 0 1 -0.1345116 51.51406 ## 201 0 1 -0.1344245 51.51408 ## 202 0 1 -0.1346607 51.51410 ## 203 0 2 -0.1348155 51.51407 ## 204 0 3 -0.1353288 51.51400 ## 205 0 3 -0.1353251 51.51376 ## 206 0 1 -0.1351935 51.51378 ## 207 0 1 -0.1351738 51.51355 ## 208 0 3 -0.1350901 51.51364 ## 209 0 1 -0.1351064 51.51345 ## 210 0 1 -0.1345174 51.51352 ## 211 0 1 -0.1343525 51.51357 ## 212 0 1 -0.1342768 51.51359 ## 213 0 2 -0.1341818 51.51362 ## 214 0 2 -0.1342435 51.51352 ## 215 0 4 -0.1344024 51.51349 ## 216 0 5 -0.1344929 51.51345 ## 217 0 2 -0.1345638 51.51344 ## 218 0 5 -0.1338792 51.51331 ## 219 0 5 -0.1337681 51.51349 ## 220 0 3 -0.1339759 51.51329 ## 221 0 3 -0.1340732 51.51326 ## 222 0 1 -0.1342086 51.51322 ## 223 0 5 -0.1342988 51.51318 ## 224 0 4 -0.1343863 51.51316 ## 225 0 4 -0.1346111 51.51309 ## 226 0 1 -0.1349730 51.51297 ## 227 0 4 -0.1350693 51.51295 ## 228 0 1 -0.1351586 51.51292 ## 229 0 3 -0.1352713 51.51289 ## 230 0 2 -0.1353472 51.51285 ## 231 0 1 -0.1356248 51.51287 ## 232 0 2 -0.1350448 51.51334 ## 233 0 1 -0.1348969 51.51337 ## 234 0 1 -0.1357609 51.51306 ## 235 0 2 -0.1358162 51.51311 ## 236 0 3 -0.1358666 51.51323 ## 237 0 1 -0.1366945 51.51341 ## 238 0 1 -0.1357570 51.51326 ## 239 0 4 -0.1363895 51.51299 ## 240 0 2 -0.1365337 51.51320 ## 241 0 2 -0.1366338 51.51313 ## 242 0 1 -0.1366662 51.51320 ## 243 0 5 -0.1364772 51.51355 ## 244 0 3 -0.1363060 51.51424 ## 245 0 2 -0.1361014 51.51428 ## 246 0 3 -0.1355020 51.51402 ## 247 0 2 -0.1354588 51.51420 ## 248 0 1 -0.1368681 51.51180 ## 249 0 1 -0.1365176 51.51149 ## 250 0 1 -0.1361566 51.51135
ggmap(m) +
geom_point(data = cholera_proj,
aes(x = coords.x1, y = coords.x2, size = Count))
Datum
The datum is the point of reference of a GIS system. To fully specify how to put data on a map, you need both the projection and the datum.
Our original cholera data didn't have a datum specified. Some digging will show that it was meant to be in the British standard, so we set that as follows.
proj4string(cholera) <- CRS("+init=epsg:27700")
## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform
Reproject
cholera_proj <- cholera %>%
spTransform(CRS("+init=epsg:4326")) %>%
as.data.frame()
ggmap(m) +
geom_point(data = cholera_proj,
aes(x = coords.x1, y = coords.x2, size = Count))
Add in pumps
proj4string(pumps) <- CRS("+init=epsg:27700")
## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform
pumps_proj <- pumps %>%
spTransform(CRS("+init=epsg:4326")) %>%
as.data.frame()
Full map code
ggmap(m) +
geom_point(data = cholera_proj,
aes(x = coords.x1, y = coords.x2, size = Count)) +
geom_point(data = pumps_proj,
aes(x = coords.x1, y = coords.x2), size = 3, color = "red", pch = 18)
John Snow's Cholera Map
Let's brush this up by using a basemap with a more similar aesthetic to the original.
m <- get_map("John Snow, London, England", zoom = 17, maptype = "toner")
Homework 7
This homework to be done in pairs (but each in your own repo). Get started.